#include "atmoslayer.h"
#include <QDataStream>
#include <QTimer>
#include <math.h>
#include <QDebug>
#include <QTime>
#include <QDir>
//qDebug() << "before segement fault.******************";
#define CONSERVE_TEST
void AtmosLayer::emitLight()
{
#ifdef CONSERVE_TEST
    int iThetaS, iPhiS;
    double intensUp, intensD;
    intensUp = intensD = 0;
    for(iThetaS = 0; iThetaS > -N_PI/2; iThetaS--){

	for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

	    intensD += (m_lightScatD[-iThetaS][iPhiS][0]) *
		    cos(iThetaS*M_PI/N_PI) *
		    sin(M_PI + iThetaS*M_PI/N_PI)*
		    pow(M_PI/N_PI,2);
	}
    }
    intensD /= M_PI; //Pi*F = Int I*u dw, F = 1/Pi Int I*u dw.

    for(iThetaS = 0; iThetaS < N_PI/2; iThetaS++){

	for(iPhiS = 0; iPhiS < 2*N_PI; iPhiS++){

	    intensUp += (m_lightScatUp[iThetaS][iPhiS][0]) *
		    cos(iThetaS*M_PI/N_PI) *
		    sin(iThetaS*M_PI/N_PI)*
		    pow(M_PI/N_PI,2);
	}
    }
    intensUp /= M_PI; //Pi*F = Int I*u dw, F = 1/Pi Int I*u dw.

    qDebug() << "Scattering irradianceUp" << intensUp << "irradianceD" << intensD << "(irradiance up+down)" << (intensUp + intensD);

#endif

    qDebug() << m_layerName << "emitting light, current caclating direction" << (m_currentDir ? "down" : "up");

    QByteArray data;
    QDataStream strm(&data, QIODevice::ReadWrite);

    if(m_currentDir == DOWN){
	//format light data to QByteArray
	for(int i = 0; i < N_PI/2; i++){
	    for(int j = 0; j < 2*N_PI; j++){
		strm << m_lightScatD[i][j][0] <<
			m_lightScatD[i][j][1] <<
			m_lightScatD[i][j][2];

		//accumulate total light as final result
		m_lightScatDTotal[i][j][0] +=
			m_lightScatD[i][j][0];
		m_lightScatDTotal[i][j][1] +=
			m_lightScatD[i][j][1];
		m_lightScatDTotal[i][j][2] +=
			m_lightScatD[i][j][2];

		//clear when sendout
		m_lightScatD[i][j][0] =
			m_lightScatD[i][j][1] =
			m_lightScatD[i][j][2] = 0;
	    }
	}

	//damped solar flux must emit first!
	if(m_solarOut > 0){
	    emit(emitSolarLight(m_solarOut,
				m_solarZen, m_solarAzi));
	    //clear when send out.
	    m_solarOut = 0;
	    qDebug() << QTime::currentTime().toString("hh:mm:ss") << m_layerName << "emitSolarLight() over.";
	}

	saveResult(DOWN);

	//scatter light emit to low layer
	//and be sure to return immediately, avoid recursion.
	//lower layer designed to start caculate after 10ms.
	emit(emitLightDown(data));
	qDebug() << QTime::currentTime().toString("hh:mm:ss") << m_layerName << "emitLightDown() over.";


    } else{
	//format light data to QByteArray
	for(int i = 0; i < N_PI/2; i++)
	    for(int j = 0; j < 2*N_PI; j++){
		strm << m_lightScatUp[i][j][0] <<
			m_lightScatUp[i][j][1] <<
			m_lightScatUp[i][j][2];

		//accumulate total light as final result
		m_lightScatUpTotal[i][j][0] +=
			m_lightScatUp[i][j][0];
		m_lightScatUpTotal[i][j][1] +=
			m_lightScatUp[i][j][1];
		m_lightScatUpTotal[i][j][2] +=
			m_lightScatUp[i][j][2];
		//debug
		if(i ==1 && j == 1)
		    qDebug() << "emittingScatUp"
			     << i << j
			     << m_lightScatUp[i][j][0] << m_lightScatUp[i][j][1] << m_lightScatUp[i][j][2] << m_lightScatUpTotal[i][j][0] << m_lightScatUpTotal[i][j][1] << m_lightScatUpTotal[i][j][2];

		//clear when sendout
		m_lightScatUp[i][j][0] =
			m_lightScatUp[i][j][1] =
			m_lightScatUp[i][j][2] = 0;
	    }

	saveResult(UP);

	//no upward solar light of course
	//scatter light emit to upper layer
	//and be sure to return immediately, avoid recursion.
	//lower layer designed to start caculate after 10ms.
	emit(emitLightUp(data));
	qDebug() << QTime::currentTime().toString("hh:mm:ss") << m_layerName << "emitLightUp() over.";
    }
}

void AtmosLayer::recvScatInDown(QByteArray lightData)
{
	qDebug() << m_layerName << "recieved downward scatter light";

	QDataStream strm(lightData);

	for(int i = 0; i < N_PI/2; i++)
		for(int j = 0; j < 2*N_PI; j++)
			strm >> m_lightInD[i][j][0] >>
				m_lightInD[i][j][1] >>
				m_lightInD[i][j][2];

	QTimer::singleShot(10, this, SLOT(downInResponse()));
}

void AtmosLayer::recvScatInUp(QByteArray lightData)
{
	qDebug() << m_layerName << "recieved upward scatter light";

	QDataStream strm(lightData);

	for(int i = 0; i < N_PI/2; i++){
		for(int j = 0; j < 2*N_PI; j++){
			strm >> m_lightInUp[i][j][0] >>
				m_lightInUp[i][j][1] >>
				m_lightInUp[i][j][2];
		}

	}

	QTimer::singleShot(10, this, SLOT(upInResponse()));
}

void AtmosLayer::recvSolar(double intensity,
			   double zenith, double azimuth)
{
	if(zenith > M_PI || zenith < M_PI_2){
		qDebug() << "Zenith must be in (Pi/2, Pi) for downward solar flux, while you have input zenith:" << zenith << ". Exiting...";
		exit(-1);
	}

	m_solarIntensity = intensity;
	m_solarZen = zenith;
	m_solarAzi = azimuth;

	qDebug() << m_layerName << "recieved solar direct light. Intensity"
		 << m_solarIntensity << "zenith" << m_solarZen
		 << "azimuth" << m_solarAzi;

	qDebug() << "irradiance loss of solar flux" <<
		    m_solarIntensity*qAbs(cos(m_solarZen))
		    *(1-exp(-1*m_tau/qAbs(cos(m_solarZen))));

}

#include <QtConcurrentRun>
void AtmosLayer::upInResponse()
{
	m_loopUpTime++;

	qDebug() << QTime::currentTime().toString("hh:mm:ss") << m_layerName << "calculating upward incident light's response for the " << m_loopUpTime << "times.";

	QFuture<void> up = QtConcurrent::run(
				this, &AtmosLayer::upInUpScat);
	QFuture<void> down = QtConcurrent::run(
				this, &AtmosLayer::upInDownScat);

	up.waitForFinished();
	down.waitForFinished();

	//when reponse caculate over, clear incident light
	for(int i = 0; i < N_PI/2; i++)
		for(int j = 0; j < 2*N_PI; j++)
			m_lightInUp[i][j][0] =
					m_lightInUp[i][j][1] =
					m_lightInUp[i][j][2] = 0;

	emitLight();

	return;
}

void AtmosLayer::downInResponse()
{
	m_loopDownTime++;

	qDebug() << QTime::currentTime().toString("hh:mm:ss") << m_layerName << "calculating downward incident light's response for the " << m_loopDownTime << "times.";

	//*************first caculate scatter-incident light response
	downScatInResponse();

	//************then caculate solar-incident light response
	if(m_solarIntensity > 0)
	    solarScatCaculate();

	//caculate downward irridiance, sendout getDownIrridiance() signal
	//must caculate here before emitLight()(downLight will be clear)
	downIrradianceCaculate();

	//************response caculate over, emit light out
	emitLight();

	return;
}
//#define SINGLE_TREAD
void AtmosLayer::downScatInResponse()
{
#ifdef SINGLE_TREAD
    downInUpScat();
    downInDownScat();
#else
    //depart to two parrallel caculating thread,
	//one for upward scatter, the other for downward scatter.
	QFuture<void> up = QtConcurrent::run(
				this, &AtmosLayer::downInUpScat);
	QFuture<void> down = QtConcurrent::run(
				this, &AtmosLayer::downInDownScat);

	up.waitForFinished();
	down.waitForFinished();
#endif
//	downInDownScat();
//	downInUpScat();

	//when reponse caculate over, clear incident light
	for(int i = 0; i < N_PI/2; i++)
		for(int j = 0; j < 2*N_PI; j++)
			m_lightInD[i][j][0] =
					m_lightInD[i][j][1] =
					m_lightInD[i][j][2] = 0;
}

void AtmosLayer::solarScatCaculate()
{
	qDebug() << "AtmosLayer::solarScatCaculate().";
	if(m_solarIntensity != 0){
		QFuture<void> solarUp = QtConcurrent::run(
					this,
					&AtmosLayer::solarUpScatCaculate);
		QFuture<void> solarDown = QtConcurrent::run(
					this,
					&AtmosLayer::solarDownScatCaculate);

		solarUp.waitForFinished();
		solarDown.waitForFinished();

		//caculate damped solar flux that will incident to next layer
		m_solarOut = m_solarIntensity *
			exp(-1*m_tau/qAbs(cos(m_solarZen)));

		m_solarIntensity = 0;

/*		qDebug() << "Symmetry:\n"
			 << 200 << m_lightScatD[45][200][0] << m_lightScatD[45][200][1] << m_lightScatD[45][200][2] << "\n"
			 << 200 << m_lightScatD[65][200][0] << m_lightScatD[65][200][1] << m_lightScatD[65][200][2];
*/
	}


}

#include <QFile>
#include <math.h>
void AtmosLayer::saveResult(bool dir)
{
    //form two file, 3-Stokes vector at all angle && degree of polarization at all angle
    QString stokesFileName, polarFileName;
    QByteArray stokesFileData, polarFileData;

    if(m_resultPath.isEmpty()){
	//format result path, mkdir
	m_resultPath.append("result/");
	m_resultPath.append(QString().sprintf("%.3fum-%.2f_%.2f-", m_wl, m_solarZen*180/M_PI, m_solarAzi*180/M_PI));
	m_resultPath.append(QDate::currentDate().toString("yyMMdd"));
    }

    QDir directory;
    if(!directory.exists(m_resultPath)){
	qDebug() << "making results' directory at " + m_resultPath;
	if(!directory.mkpath(m_resultPath)){
	    qDebug() << "!!!!!!!!!!!!!!!!!!!!!!!!!\nCannot make results' directory at " + m_resultPath + "\n!!!!!!!!!!!!!!!!!!!!!!!\n";
	    exit(-2);
	};
    }

    stokesFileName = polarFileName = m_resultPath;

    if(dir == DOWN){
	if(!(m_saveMode & SAVE_DOWN))
	    return;

	stokesFileName.append(QString().sprintf("/%.2fm_down_%d.txt",
						m_hBottom, m_loopDownTime));
	polarFileName.append(QString().sprintf(
				 "/%.2f_down_%d_polar_degree.txt",
				 m_hBottom, m_loopDownTime));

	//form first azimuth index line.
	QString stokesFirstLine("Zenith\\Azimuth\t");
	QString polarFirstLine("Zenith\\Azimuth\t");
	int aziIndex;
	double azimuth;
	for(aziIndex = 0, azimuth=0;
	    aziIndex < 2*N_PI; aziIndex++,
	    azimuth += ANGLE_STEP){

	    stokesFirstLine.append(QString().sprintf(
				       "I_a_%.2f\tQ_a_%.2f\tU_a_%.2f\t",
				       azimuth*180/M_PI, azimuth*180/M_PI, azimuth*180/M_PI));
	    polarFirstLine.append(QString().sprintf(
				      "P_a_%.2f\t",
				      azimuth*180/M_PI));
	}
	stokesFirstLine.append("\n");
	stokesFileData.append(stokesFirstLine);

	polarFirstLine.append("\n");
	polarFileData.append(polarFirstLine);

	QString stokesZenithLine, polarZenithLine;
	double zenith = M_PI - ZENITH_INIT; //[Pi, Pi/2]
//	double zenith = M_PI; //[Pi, Pi/2]
	for(int i = 0; i < N_ZENITH; i++, zenith -= ANGLE_STEP){

	    //first word is zenith angle index
	    stokesZenithLine.append(QString().sprintf("Z%.2f\t",
						      zenith*180/M_PI));
	    polarZenithLine.append(QString().sprintf("Z%.2f\t",
						     zenith*180/M_PI));

	    for(int j = 0; j < 2*N_PI; j++){

		stokesZenithLine.append(QString().sprintf(
					    "%f\t%f\t%f\t",
					    m_lightScatDTotal[i][j][0],
					    m_lightScatDTotal[i][j][1],
					    m_lightScatDTotal[i][j][2]));
//		if(i == 0)
//		    qDebug() << i << j
//			     << m_lightScatDTotal[i][j][0]
//			     << m_lightScatDTotal[i][j][1]
//			     << m_lightScatDTotal[i][j][2];
		/*
    //degree of polarization = sqrt(Q^2 + U^2) / I
    polarZenithLine.append(QString().sprintf(
    "%f ",
    (sqrt(pow(m_lightScatDTotal[i][j][1], 2) +
     pow(m_lightScatDTotal[i][j][2], 2))) /
     m_lightScatDTotal[i][j][0]));
*/
		polarZenithLine.append(QString().sprintf(
					   "%f\t",
(sqrt(pow(m_lightScatDTotal[i][j][Q],2)+pow(m_lightScatDTotal[i][j][U],2))/m_lightScatDTotal[i][j][I])));

	    }
	    //current zenith's azimuth loop over
	    stokesZenithLine.append("\n");
	    stokesFileData.append(stokesZenithLine);//write to file cache
	    stokesZenithLine.clear(); //clear for next zenith angle
	    polarZenithLine.append("\n");
	    polarFileData.append(polarZenithLine);

	    polarZenithLine.clear();
	}

    } else{

	if(!(m_saveMode & SAVE_UP))
		    return;

	stokesFileName.append(QString().sprintf("/%.2fm_up_%d.txt",
						m_hTop, m_loopUpTime));
	polarFileName.append(QString().sprintf(
				 "/%.2f_up_%d_polar_degree.txt",
				 m_hTop, m_loopUpTime));

	//form first azimuth index line.
	QString stokesFirstLine("Zenith\\Azimuth\t");
	QString polarFirstLine("Zenith\\Azimuth\t");
	int aziIndex;
	double azimuth;
	for(aziIndex = 0, azimuth=0;
	    aziIndex < 2*N_PI; aziIndex++,
	    azimuth += ANGLE_STEP){
	    stokesFirstLine.append(QString().sprintf(
				       "I_a_%.2f\tQ_a_%.2f\tU_a_%.2f\t",
				       azimuth*180/M_PI, azimuth*180/M_PI, azimuth*180/M_PI));
	    polarFirstLine.append(QString().sprintf(
				      "P_a_%.2f\t", azimuth*180/M_PI));

	    /*			stokesFirstLine.append(QString().sprintf("A%f   ",
	azimuth*180/M_PI));
   polarFirstLine.append(QString().sprintf(
	  "A%f ",
	azimuth*180/M_PI));
*/		}
	stokesFirstLine.append("\n");
	stokesFileData.append(stokesFirstLine);

	polarFirstLine.append("\n");
	polarFileData.append(polarFirstLine);

	QString stokesZenithLine, polarZenithLine;
	double zenith = ZENITH_INIT; //[0, Pi/2]
//	double zenith = 0; //[0, Pi/2]
	for(int i = 0; i < N_ZENITH; i++, zenith += ANGLE_STEP){

	    //first word is zenith angle index
	    stokesZenithLine.append(QString().sprintf("Z%.2f\t",
						      zenith*180/M_PI));
	    polarZenithLine.append(QString().sprintf("Z%.2f\t",
						     zenith*180/M_PI));

	    for(int j = 0; j < 2*N_PI; j++){

		stokesZenithLine.append(QString().sprintf(
					    "%f\t%f\t%f\t",
					    m_lightScatUpTotal[i][j][0],
					    m_lightScatUpTotal[i][j][1],
					    m_lightScatUpTotal[i][j][2]));

		//degree of polarization = sqrt(Q^2 + U^2) / I
		polarZenithLine.append(QString().sprintf(
					   "%f\t",
(sqrt(pow(m_lightScatUpTotal[i][j][Q],2)+pow(m_lightScatUpTotal[i][j][U],2))/m_lightScatUpTotal[i][j][I])));

	    }
	    //current zenith's azimuth loop over
	    stokesZenithLine.append("\n");
	    stokesFileData.append(stokesZenithLine);//write to file cache
	    stokesZenithLine.clear(); //clear for next zenith angle
	    polarZenithLine.append("\n");
	    polarFileData.append(polarZenithLine);
	    polarZenithLine.clear();
	}

    }

    QFile file(stokesFileName);
    file.open(QIODevice::WriteOnly);
    file.write(stokesFileData);
    qDebug() << "Result saved as " + file.fileName();
    file.close();

    file.setFileName(polarFileName);
    file.open(QIODevice::WriteOnly);
    file.write(polarFileData);
    file.close();
}
